Data Preparation

# smacof analysis of 1988 wine data
# Lew Harvey
# 3 February 2016

library("smacof")
library("rgl")

Read in the data and wine properties

fn <- "wine_88.csv"
df <- read.csv(fn, header = TRUE)
df$w10 <- as.numeric(df$w10)

# read in wine properties
df.prop <- read.csv("wine_properties.csv", header = TRUE)
df.prop$wine <- as.character(df.prop$wine)
df.prop$color <- as.character(df.prop$color)

The data in the input file looks like this for the first subject:

head(df, 10)

##    w01 w02 w03 w04 w05 w06 w07 w08 w09 w10
## 1    0  NA  NA  NA  NA  NA  NA  NA  NA  NA
## 2    2   0  NA  NA  NA  NA  NA  NA  NA  NA
## 3    2   2   0  NA  NA  NA  NA  NA  NA  NA
## 4    4   4   2   0  NA  NA  NA  NA  NA  NA
## 5    2   2   6   2   0  NA  NA  NA  NA  NA
## 6    5   4   6   4   7   0  NA  NA  NA  NA
## 7    3   5   7   7   5   3   0  NA  NA  NA
## 8    5   6   7   6   5   3   7   0  NA  NA
## 9    3   5   6   6   4   2   6   5   0  NA
## 10   2   4   3   3   7   5   6   6   2   0

Convert the data into a list of matrices, one for each subject. This is a bit more complicated than it should be, but the original data are in an SPSS file that requires some massaging to get the data into a form for R

# loop through the data frame to
# create separate matrices for the
# data sets of each subject

nstim <- ncol(df)
nsets <- nrow(df) / ncol(df)
mlist <- list(NA)
mm <- list(NA)
winename <- colnames(df)

for (i in 1:nsets){
    mm[i] <- paste("m", i, sep = "")
    i1 <- ((i - 1) * nstim) + 1
    i2 <- i1 + nstim - 1
    mlist[i] <- list(mm = as.matrix(df[i1:i2,]))
    row.names(mlist[[i]]) <- winename
    mlist[[i]] <- as.dist(mlist[[i]])
}
names(mlist) <- unlist(mm)

Now the matrices are in a list, one for each subject. Here are the first subject’s data as a distance matrix:

mlist[1]

## $m1
##     w01 w02 w03 w04 w05 w06 w07 w08 w09
## w02   2                                
## w03   2   2                            
## w04   4   4   2                        
## w05   2   2   6   2                    
## w06   5   4   6   4   7                
## w07   3   5   7   7   5   3            
## w08   5   6   7   6   5   3   7        
## w09   3   5   6   6   4   2   6   5    
## w10   2   4   3   3   7   5   6   6   2

2D INDSCAL Solution

wine2D <- smacofIndDiff(mlist,
                      ndim = 2,
                      type = "ordinal",
                      constraint = "indscal",
                      itmax = 1000)
wine2D

## 
## Call: smacofIndDiff(delta = mlist, ndim = 2, type = "ordinal", constraint = "indscal", 
##     itmax = 1000)
## 
## Model: Three-way SMACOF 
## Number of objects: 10 
## Stress-1 value: 0.2201172 
## Number of iterations: 440

Plot 2D Solution and Fits

par(mfcol = c(2, 2))
par(pty = "s")
plot(wine2D, type = "p",
     plot.type = "confplot",
     label.conf = c(TRUE, 1, 1),
     xlim = c(-0.9, 0.9),
     ylim = c(-0.9, 0.9))
abline(h = 0, col = "gray")
abline(v = 0, col = "gray")
par(pty = "m")

plot(wine2D, plot.type = "Shepard") 
plot(wine2D, plot.type = "stressplot", ylim = c(0, 25))
plot(wine2D, plot.type = "resplot")

par(mfcol = c(1, 1))

configuration plot without color

par(pty = "s")
plot(wine2D, type = "p", pch = 19,
     plot.type = "confplot",
     label.conf = c(TRUE, 1, 1),
     xlim = c(-0.9, 0.9),
     ylim = c(-0.9, 0.9))
abline(h = 0, col = "gray")
abline(v = 0, col = "gray")

par(pty = "m")

configuration plot with color

par(pty = "s")
plot(wine2D, type = "p", pch = 19,
     col = df.prop$color,
     plot.type = "confplot",
     label.conf = c(TRUE, 1, 1),
     xlim = c(-0.9, 0.9),
     ylim = c(-0.9, 0.9))
abline(h = 0, col = "gray")
abline(v = 0, col = "gray")

par(pty = "m")

3D solution

wine3D <- smacofIndDiff(mlist,
                        ndim = 3,
                        type = "ordinal",
                        constraint = "indscal",
                        verbose = FALSE,
                        itmax = 5000)
wine3D

## 
## Call: smacofIndDiff(delta = mlist, ndim = 3, type = "ordinal", constraint = "indscal", 
##     verbose = FALSE, itmax = 5000)
## 
## Model: Three-way SMACOF 
## Number of objects: 10 
## Stress-1 value: 0.155876 
## Number of iterations: 1688

Plot the 3D solution without color

This plot cannot be printed since it is 3D and can be rotated by the user

plot3d(wine3D$gspace, type = "s",
       expand = 1.3,
       xlab = "X",
       ylab = "Y",
       zlab = "Z")
text3d(wine3D$gspace,
       text = df.prop$wine,
       adj = c(1.5, 1.5))

Plot the 3D solution with color